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Abstract — Computational techniques are proposed for a numerical solution of an inverse problem that arises in a study of a 
microbial depolymerization processes. The Newton-Raphson method in conjunction with the Newton ’s method is applied to 
a time factor that involves three parameters. The Newton-Raphson method reduces a system of three equations for three 
unknowns to a single variable equation, which is solved with the Newton ’s method. Those techniques were applied to a 
microbial depolymerization process of polyethylene glycol. Introduction of experimental results into analysis leads to 
simulation of a microbial depolymerization problem. 
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I. Introduction 

Petroleum-based polymers have been produced since the early part of the twentieth century. A large amount of those 
xenobiotic compounds have been accumulated on the surface of the earth, and they are now sources for carbon dioxide 
emission. Although those macromolecules had been nonexistent until they were invented, some of microorganisms are able 
to utilize the carbon sources. Mechanisms of microbial depolymerization processes must be clarified for appropriate 
assessment of the carbon cycle. Microbial depolymerization processes are categorized into exogenous type processes and 
endogenous type processes. Monomer units are liberated one by one from terminals of molecules in an exogenous type 
depolymerization process. Polymers degraded in exogenous type depolymerization processes include polyethylene (PE) and 
polyethylene glycol (PEG). Unlike exogenous type depolymerization processes, primary factor of an endogenous type 
depolymerization process is arbitrarily scission. Polymers degraded in endogenous type depolymerization processes include 
polyvinyl alcohol (PVA) and polylactic acid (PLA). 

A PE biodegradation process involves two primary factors, gradual weight reduction due to the p -oxidation and direct 
absorption by cells. A PE molecule liberates a two carbon unit from its terminal in one cycle of p -oxidation and reduces in 
size undergoing successive p -oxidation processes until it becomes small enough to be absorbed directly into cells. A 
mathematical model based on those factors was proposed for simulation of PE biodegradation processes [1 - 4], 
Polyethylene glycol (PEG) is another polymer depolymerized in exogenous type depolymerization processes. The initial 
step of liberation from a PEG molecule is oxidation of a terminal unit to produce an aldehyde, which is further oxidized to a 
monocarboxylic acid. The oxidation is followed by the cleavage of the ether bond, and the molecule reduces by one glycol 
unit [27]. Numerical techniques developed for PE biodegradation were applied for simulation of an exogenous type 
depolymerization process of PEG [5]. Temporal factors of degradation rates were taken into consideration [6-8], 

A mathematical model of endogenous type depolymerization processes was proposed for an enzymatic degradation process 
of PVA [9, 10]. The model was applied to an enzymatic hydrolysis of polylactic acid (PLA), and degradabilities of PVA and 
PLA were compared [11], Temporal dependence of the degradation rate was taken into consideration for the 
depolymerization process of PLA [14]. The model proposed for endogenous type depolymerization processes was 
reformulated for exogenous type depolymerization processes of PEG [13] and PE [14]. Numerical techniques developed for 
PE biodegradation were applied to an exogenous depolymerization process of PEG [15], Temporal dependence of 
degradability was taken into consideration for depolymerization processes of PEG [16 - 23, 25, 26]. 

The study on microbial depolymerization processes of PEG was continued. In previous studies, weight distributions of PEG 
before and after cultivation of microbial consortium El were introduced into analysis, and techniques to evaluate parameters 
for a time factor of degradation rate were proposed. In this study, those techniques were improved. The Newton-Raphson 
method reduces a system of three equations with three unknowns to a single equation of one variable equation which is 
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solved numerically by the Newton's method. Our techniques are illustrated with results obtained from application to weight 
distributions before and after cultivation of Microbial consortium El for two days, four days, and seven days 


II. Description of problems in analysis of exogenous type depolymerization processes 

Molecules of a single polymer all have the same chemical structure, while their molecular weight varies from one molecule 
to another. In order to model a microbial depolymerization process, let t and M be the time and the molecular weight, 
respectively, and let w(t, M) be the weight distribution with respect to the molecular weight M at time t , that is, the total 
mass of polymer molecules for A< M < B is the integral 

j B w(t,M)dM. (1) 


Let a(t) be the microbial population at time t . The following system of equations was proposed in previous studies ( cf. 
[25]) under the assumption that the liberated amount distributes exponentially [24]. 


™ = cr{t\- X(M)w + c{M)^ M X(K)d{K)w{t,K) dK , 


( 2 ) 


da 

dt 


= k\ 


1 -h 



(3) 


where 


c{M)=Me pM , d{K ) 


- P K 


&e-~y 


Function V 


(0 


is the residual in the culture medium at time t , that is. 


v(f) = | w(t,M)dM, 

which is approximated as the integral (1) for appropriate values of A and B . The microbial consortium El is incubated in a 
culture medium where PEG 6000 was the sole carbon source. The weight distributions are obtained before and after 
cultivation for two days, four days, and seven days. Table 1 shows approximate values of the residual PEG y(o), v(2), v(4), 

and v’(7 ) for A = 10 32 , B = 10 42 . The system of equations (2), (3) is associated with the initial condition 

w(0 ,M)=f 0 (M ), a(0) = <r 0 , (4) 

Where / 0 (m) and a 0 are the initial weight distribution and the initial microbial population, respectively. The system of the 
equations (2), (3) and the initial condition (4) form an initial value problem. Given the initial weight distribution f 0 (M ), the 
initial value problem (2), (3), (4) is solvable provided the molecular factor A(m) and the initial microbial population cr 0 , and 
the values of the parameters k and h are given. In what follows, techniques for determination of those unknowns are 
described. Suppose that weight distributions at t = t 0 , t = t t , t = t., , and t — t 3 are expressed with prescribed functions 
/ 0 (m), /,(m), and / 3 (m), respectively. Here t 0 = 0 and / 0 (M) is the initial weight distribution. A problem 

posed in previous studies is to determine the values of the parameters cr 0 , k and h , which satisfy 

w(f 1 ,M)=f 1 (M), w(t 2 ,M)=f 2 (M), w(f 3 ,M) = f 3 (M ). (5) 
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Equation (2) involves a(t) while equation (3) involves v(f), the integral of w(t,M ) with respect to M . In order to remove 
the mutual dependence, consider the change of independent variables from f to r by 

r = j\r(.s)t/s. (6) 

Note that 

at 


holds. Let 

W(t,M)= w(t,M), S(r)=cr(t). 


In view of the expression (7), equations (2) and (3) are converted to 


= -Z(m )W + c(M )£ Z(K)d(K)w(T,K)dK, 


( 8 ) 


dS . , S 
— = 1 -n-rr, 
dr v(t) 


(9) 


Respectively. Note that the equation (8) now involves only one of the unknown variables W . which makes the inverse 
problem for the molecular factor X{M ) solvable, while the equation (2) involves two unknown variables w and a . 

Given the molecular factor of the degradation rate A(M ) and a weight distribution fAm) for r = 7’ , equation (8) and the 
initial condition 

w(t v m)=f 1 (m ) (10) 

form an initial value problem. Given an additional weight distribution fAm) for t = T . equation (8), the initial condition 

(10) , and the additional condition 

w(t v m)=f 1 (m) (11) 

form an inverse problem for X{m), for which the solution of the initial value problem (8), (10) also satisfies the condition 

(11) . Techniques to solve the inverse problem were developed in the previous studies. In particular, the molecular factor 
X (M ) was evaluated for T t = 0 and T 2 = 2 . for which F t (M ) and F n ( M ) were the weight distributions of PEG after 
cultivation of the microbial consortium El for two days and four days, respectively. Once A(lVl ) was obtained, the initial 
value problem consisting of the equation (8) and the initial condition 


W(O l ,M)=f 0 (M) 


( 12 ) 
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FIGURE 1: Transition of residual PEG for 0<T<10 based on the solution of the initial value problem (8), (12). The 
initial value problem (8), (12) was solved numerically, and curve (X, y(T)) was generated | y M = ^( r ) x | qq 


V(0) 


y(r)= \ B W(T,M)dM, A = 10 32 and B = 10 42 . 

JA 


Table 1 

Values of t and r , and residual PEG before and after cultivation of microbial consortium El for two days, four days, 
and seven days ( A = 10 3 2 and B = 10 4 2 ). The table shows the values of t , approximate values of z , and approximate 

values of the residual PEG (%) based on the conditions (5). 


t 

z 

Residual PEG 

0.00000 

0.00000 

100.00000 

2.00000 

0.51745 

72.39508 

4.00000 

2.01635 

26.59018 

7.00000 

7.21585 

1.89424 


was solved, and the residual PEG V(r) as a function of r , 


V(z)=^W(z,M)dM (13) 

was obtained (A = 10 32 and B = 10 42 ). Figure 1 shows the curve (r, V / (r)) for 0 < r < 1 0 . Values of r corresponding to 
values of t were obtained by solving the equations V(r) = i'(?) for t = 0, 2, 4, and 7 (Table 1). Table 1 shows that 
t = 2.0, t = 4.0, and t = 7.0 correspond to z ~ 0.51745, z ~ 2.01635, and z ~ 7.21585, respectively. 


III. Numerical techniques for evaluation of parameters of time factor 

The solution S of equation (9) with the initial value a 0 depends not only on r but also on cr 0 , k and h . Denote by 
S — S(z,(T 0 ,k,li) the solution of (9) which satisfies the initial condition S{0,a n ,k,li)— o Q . Let 
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u(T,CT 0 ,k,h)= f 


S{q,a 0 ,k,h) 


dq. 


In view of the expression (7), t — u(r,cr 0 ,k,h). Suppose that / , l\ , and t { correspond to r, , r 2 , and r 3 according to the 
change of variables (6), that is. 


t, =u(T i ,<J 0 ,k,h) 


0 = 1, 2, 3). 


(14) 


Let g i (cr Q ,k,h)= u(j j ,<T 0 ,k,h) — t i , and consider the system of equations for the unknowns cr 0 , k and h obtained from 
the conditions (14), 


g i (cr 0 ,k,h) = 0 0 = 1, 2, 3). 

In order to solve the system of equations (15), consider the first two equations 

g l {a n ,k,h) = 0, g 2 {a Q ,k,h) = 0 


(15) 


(16) 


for a fixed non-negative value of h . System of equations (15) can be solved for cr 0 and k numerically by the Newton- 
Raphson method to generate functions <j)(h) and i//(li) such that a solution of the system is given by cr 0 = </>(h), k — f/x(/t) , 
that is, a sequence of approximate solutions (cr 0 (| ,k 0 ), (cr 0 , k t ], (cr () 2 , Ar 2 ), is generated by the recurrence formula 


cr, 


0,m+l ^0,» 


k - k - 

K m + 1 — *,« 


f 

dg 2 

og 2 

dg' 

-1 


dg\ 

o 

b 

CO 

dk 

© 

b 

<o 

dk j 


t dk gx 

dk 


( dgi 

dg 2 

Sg 2 


-y 

© 

b 

<o 

dk 

© 

b 

<o 

dk j 

V 


Sg 2 , dg 1 

T gl +~ g 2 

8u 0 da 0 J 


m = 0,1,2, ..., 


where the functions g t and g, and their partial derivatives are evaluated at a 0 = a 0n and k =k n - Note that 


dg f du 

dg l du 

dgj du 

da 0 da 0 

’ ~dk~dk 

T=T: 

r=r ’ dh dh 


while the partial derivatives du / da n , du / dk , and on / dh are the solutions of the initial value problems 


d_ 

dr 


' du ^ 

1 dS 

du 

v 5cr o J 

[5(r)] 2 da 0 ’ 

da 0 

I f du 

) 1 dS 

du 


T 

) “ [s(r)] 2 dk ’ 

dk 

T=0 

f 

1 dS 

du 


\Wij 

[s{t)P a* • 

dh 

r=0 


= 0, 


= 0, 


= 0 . 


Partial derivatives dS / da n , dS / dk , and dS / dh are the solutions of the initial value problems 


= 1 , 


d 

( dS ^ 

kh 

r ds ^ 

dS 

dr 


V(r)l 


dc 7 0 


(17) 


(18) 


(19) 


( 20 ) 


( 21 ) 


( 22 ) 
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d_ 

dr 


dr 


f dS 'j _ hS( r) kh f dS 'i 

V3A: J V(r) V(r) v 3fcJ’ 

V(r)L W dh) 


as 


= 0, 


as 

5/t 


= o. 


Those initial value problems are solved for the partial derivatives which appear in the recurrence formulas (17) [26]. 
Note that the system (16) leads to 


</>'{ h ) = 


y/'{h) = 


dgi dg 2 d § 2 d 8 i 

v 5cr () dk 8cr 0 dk y 

dgi dg 2 dg 2 dg x 

v a<7 0 a^ a<T (l a^ 


dg 2 ggi fjgi % 2 
<3& dh dh dk 

dg 2 dg l dg l dg 2 


yda 0 dh da 0 dh 


for cr 0 =0(h) and k=if/(h). Let 


Z(h) = g :t (d(h),i//(h),h). 

Once the solution <j 0 = </>(h), k — !//(/;] of the system (16) is obtained, the Newton's method can be applied to 

Z(h)= 0, 


(23) 


(24) 


(25) 


(26) 


(27) 


that is, a sequence of approximate solutions h 0 ,h t ,h n can be generated by 


Note that 


K + i = K - 


Z(h„) 

Z'ihJ 


n = 0, 1, 2, 


' ' rs _ ' rsi. 1 rsi 


da. 


dk 


dh 


(28) 


for cr n = (f){h ) and k = y/{h ). Derivatives (j)' and y/' are given by the expressions (25), while the partial derivatives of g ^ , 
g , , and g 3 with respect to a 0 , k and h are given by expressions (18) - (24). 


IV. Numerical results from analysis of inverse problem for time factor 

Let V(t) be the residual PEG obtained from a numerical solution of the initial value problem (8), (12), that is, V(t) is the 
function defined by the expression (13) for A = 10 3 2 and B = 10 42 . The system of equations (15) was solved numerically for 
cr 0 , k and h with the techniques described in the previous section. Table 2 shows the results from application of the 

Newton-Raphson method to system (16) and the Newton's method to equation (27). Figure 2 shows the transition of the 
residual PEG obtained from a numerical solution of the initial value problem (8), (9), (12) with 5 ( 0 , a 0 , k, h) = a 0 ■ Figure 2 

also shows the transition of residual PEG obtained from a numerical solution of the initial value problem (2), (3), (4) for the 
values of the parameters a Q , k and h in the last row of table (2). The initial value problem (2), (3), (4) was solved 

numerically using values of parameters a Q , k and h in the last row of table (2). Figure 3 shows the transition of the residual 
PEG obtained from a numerical solution of (2), (3), (4). 
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Table 2 

Change of parameter values and iteration counts of the Newton -Raphson method and the Newton's method. The 
table shows that /; ] - h n | < i o -12 for n = 5 in the sequence generated by the Newton's method (28). It also shows that 

V( CT o.„, + 1 ~ a o, m ) 2 + {k m+ i-k m T ^ 10 12 for m = 16 at each step of the Newton's method. Entries in the columns 
cr 0 n and k n are the last values of the Newton-Raphson iterations at the n th step of the Newton's method. The initial 
value of each application of the Newton-Raphson method was fixed ( cr 0 0 = 0.2 , k 0 = 0.5 ). 


n 

K + , 

Count 

^0 ,n 

k n 


0 

12.91011526876690 

16 

0.1 4381 1 26167661 

0.54194442446577 

-0.20163026883843 

1 

12.75353240846410 

16 

0.1 43432045941 1 6 

0.54498679389487 

0.01215719681922 

2 

1 2.7531 631 5092790 

16 

0.14345256753083 

0.54482207592820 

0.00002849758473 

3 

1 2.7531 6339202850 

16 

0.14345261590947 

0.54482168762418 

-0.00000001860676 

4 

12.75316339186900 

16 

0.1 4345261 587788 

0.54482168787771 

0.00000000001231 

5 

12.75316339186950 

16 

0.1 4345261 587790 

0.54482168787755 

-0.00000000000004 
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FIGURE 2: Transition of residual PEG for seven days based on the initial value problem (8), (9), (12). The curve 


(u(t, a 0 ,k , /?), y{r)) is shown 


^)-yi xl0 ° 


Function u(r , <j 0 ,k , h ) was generated from a numerical solution of the 


initial value problem (8), (9), (12) with the initial condition ,S'(0, a 0 ,k , h ) = a 0 for the values of the parameters cr 0 , k 

and h in the last raw of table 2. 
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FIGURE 3: Transition of residual PEG for seven days based on the initial value problem (2), (3), (4). The curve 
( t,z(t )) is shown ( _M_ V v) ;; iqq | • Function v(l ) was generated from a numerical solution of the initial value 

r v v (°) J 

problem (2), (3), (4) for the values of the parameters cr 0 , k and h in the last raw of table 2. 

Figure 1 shows that the residual PEG V (r) may well be approximated by an exponential function. Suppose that V (r) is a 
solution of 


(29) 


(30) 


dV 

so that V(t) — V(o)e~ f,T . Suppose that the value of u is determined by the conditions 

V(r,) = V 1 , V(r 1 ) = V 2 , 

where Vj and V 1 are the residual PEG at r — and v — T 1 , respectively, so that 

logE.-logF, 

r 2 -r. 

Figure 4 shows the curve (r,e~ MT x lOO). The figure also shows the curve (r, y(r)) (Figure 1). Parameters <j Q . k and /? 
were evaluated for = V' () e' /,r with the techniques described in the previous section, where V 0 is the residual PEG before 
cultivation, that is. 


y o =jX(0 ,M)dM (a = 10 32 , S = 10 42 ). 
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A = 10 


and B = 1 O 4 2 . 


Table 3 


Change of parameter values and iteration counts of the Newton-Raphson method and the Newton's method for 
V(r)=V 0 e~ f ‘ I • The table shows that \h n+i -h n \<10~ 12 for n= 5 in the sequence generated by the Newton's method (28). 

It also shows that V( a o,,„ + i -CT 0 J 2 + (k m+l — k m ) 2 <10 12 for m = 16 at each step of the Newton's method. . 
Entries in the columns <j 0 n and k n are the last values of the Newton-Raphson iterations at the n th step of the 
Newton's method. The initial value of each application of the Newton-Raphson method was fixed 

(o- 00 =0.2 ,k 0 =0.5). 


n 

K + 1 

Count 

^O.n 

k n 

83(<*0j» k n>K) 

0 

8.58671318397139 

16 

0.14379681030420 

0.542081 1 1 099523 

0.18837377402843 

1 

8.51606276891822 

16 

0.14398144583988 

0.54059785942797 

0.00856257414694 

2 

8.51591999712240 

16 

0.14399064713205 

0.54052396031086 

0.00001721986678 

3 

8.51592005463291 

16 

0.14399066572345 

0.54052381099793 

-0.00000000693634 

4 

8.51592005460980 

16 

0.1 4399066571 597 

0.54052381105807 

0.00000000000279 

5 

8.51592005460968 

16 

0.1 4399066571 597 

0.54052381105804 

0.00000000000001 


Figure 5 shows the transition of the residual PEG obtained from a numerical solution of the initial value problem (8), (9), 
(12) with s(0,a 0 ,k,h)= a 0 and v(r) = V 0 e~ f ‘ T for the values of the parameters cr 0 , k and h in the last raw of table 3. 
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FIGURE 5: Transition of residual PEG for seven days based on the initial value problem (8), (9), (12). The curve 


(u(r, (T a , k, li ), y(r)j is shown 


y(r) = ^|]xl00, V(r) = V 0 e~ 


Function u(z, a 0 ,k, h) was generated from a 


numerical solution of the initial value problem (8), (9), (12) with the initial condition .S'(6, <x 0 , k, h) = cr 0 for the values 

of the parameters cr 0 , k and h in the last raw of table 3. 


V. Conclusion 

Application of the Newton-Raphson method to the system (16) was demonstrated in a previous study, while the bisection 
method was applied to the equation (27) [26]. The techniques were tested by application to model problem in which V(z) 

was a constant function. The techniques were also applied to an exponential function V(r) =e~ aT+b where the least square 
approximation was applied to V(r) based on weight distributions before and after cultivation of microbial consortium El for 
one day, three days, five days, and nine days, and the values of the parameters <x 0 ~ 0.14345, k ~ 1.5690 and h « 9.3784 
were obtained [26]. 

In this study, the application of the Newton's method to the equation (27) was illustrated. The techniques were applied to the 
residual PEG V(r) based on the weight distributions before and after cultivation of microbial consortium El for two days, 

four days, and seven days. The values of the parameters a 0 ~ 0.14345, k « 0.54482 and h « 12.75316 (Table 2) were 
obtained. Values of a 0 between numerical results in the previous study and this study are almost identical. However there 
are notable differences in values of k and li . 

Equations (16) and (27) were also solved for the solution of the equation (29) with the initial value y o , and the values of the 
parameters cr 0 ~ 0.14399, k » 0.54052 and h « 8.5 1 592 (Table 3) were obtained. Difference in values of cr 0 and k 

between shown in Table 2 and 3 is negligible. However there is a notable difference between the values of li shown in 
Tables 2 and 3. Nevertheless, transitions of residual PEG shown in Figures 3 and 5 are indistinguishable, and equation (8) 
may well be replaced with equation (29). In practical applications, there are cases where amounts of residual polymer is 
available at discrete points in time instead of weight distributions. Techniques developed in this study will be applicable to 
those cases. 
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